Mitogenomes revealed the history of bison colonization of Northern Plains after the Last Glacial Maximum

American bison demonstrated differential patterns of extinction, survival, and expansion since the terminal Pleistocene. We determined population dynamics of the Northern Great Plains bison using 40 mitochondrial genomes from radiocarbon dated remains with the age ranging from 12,226 to 167 calibrated years before present. Population dynamics correlated with environmental and anthropogenic factors and was characterized by three primary periods: terminal Pleistocene population growth starting 14,000 years ago, mid Holocene demographic stability between 6700 and 2700 years ago, and late Holocene population decline in the last 2700 years. Most diversification of mtDNA haplotypes occurred in the early Holocene when bison colonized new territories opened by retreating ice sheets. Holocene mtDNA lineages were not found in modern bison and lacked association with archaeological sites and morphological forms.

the nineteenth century anthropogenic bottleneck 1 . That decline was linked to environmental changes, but the specific alterations preceding the dramatic population decline were not fully described. Another study, based on a subset of steppe bison mtDNA sequences, demonstrated a conflicting scenario, with population growth from 50,000 to 35,000 ya, followed by climate-caused decline in N e that reached its five-fold minimum 10,000-11,000 ya. Smaller population growth was detected after 10,000 ya, reaching a population maximum 5000 ya 10 .
More sensitive demographic analysis of the same dataset 1 based on the Bayesian skyline method revealed that between 70,000 and 25,000 ya the population size was high and stable followed by a later population decline that began 25,000 ya in the beginning of the LGM and reached a N e minimum approximately 10,000 ya, at a time when many megafaunal taxa went extinct throughout North America 16 . The latter timepoint (25,000 ya) for population decline was further supported with recent data indicating that B. priscus was abundant on Alaska's North Slope 45,000-35,000 ya 18 . That study further supports that population growth then occurred from the early Holocene approximately 10,000 ya through the mid Holocene to as recently as 5000 ya, when the population reached its Holocene N e maximum. This maximum was substantially less than the population peak before the LGM and was followed by a small population decline. The reasons for bison population reduction in the late Holocene were not clearly explained except that humans could have played a role in the decline 16 .
However, conclusions based on bison in Beringia should not be extended to contemporaneous bison populations elsewhere in North America, because of ecological plasticity and genetic discontinuity of the species across continental scales 1,4,19 . Further, both studies 1,10 were based on short mtDNA sequences, with limited samples from the Holocene analyzed in context of late Pleistocene Beringian specimens from Siberia and Alaska.
Based on paleontological and archaeological records, bison population sizes in the Holocene varied across their distributional range including the Great Plains 20 . However, exact demography as the dynamics of effective population size in local populations in the past remains unknown but can be reasonably documented based on ancient mitogenome sequences available from throughout their spatial range. Bayesian demographic analyses are useful in identifying timing of changes of bison populations which can then be associated with other drivers of environmental events. To date, only 31 complete or near complete mitochondrial genome sequences from steppe bison (B. priscus) were described from eastern Beringia and the migration corridor between Cordilleran and Laurentide ice sheets 4,21,22 . Only four of them came from Holocene specimens, including one from Charlie Lake Cave (British Columbia; KX269118), two from Clover Bar Sand and Gravel (Alberta; KX269120, KX269143) and one from Whitehorse (Yukon, MF134653). No mitogenome sequences, however, are available from the expansive territory of the Holocene range, including most of the Great Plains.
The Great Plains represents the largest biome in North America with diverse climatic factors and vegetation across prairie systems 13 . To minimize effects of differential ecology on observed genetic variation, we focused on Holocene bison specimens with determined radiocarbon dates 19 from a relatively small area of the Northern Great Plains (NGP) 23 . In this study, we tested the hypothesis whether Northern Plains bison demonstrated distinct population dynamics from Beringian bison 1,10,16 , using sequences representing a global bison population containing late Pleistocene, Holocene and modern mitogenomes. To test this hypothesis, we collected and reconstructed complete mitochondrial genomes from different Holocene periods to evaluate the mostly unknown genetic variation of Holocene bison. We also tested the time of origins of mtDNA haplotypes, their survival through the Holocene, and their relationship with B. priscus and modern mtDNA lineages that survived the anthropogenic bottleneck in the nineteenth century. Finally, we estimated N e for bison roaming the NGP throughout the Holocene, including pre-extirpation N e and compared the N e values with N e in the global population.

Results
Phylogenetic analysis of mitochondrial genomes. We carried out the target enrichment of DNA libraries prepared from 47 bison bone and tooth specimens for bison mitochondrial genome sequences followed by high throughput sequencing. Forty specimens from 18 archaeological sites ( Fig. 1; Supplementary Figs. S1, S2) yielded complete mitochondrial genome sequences (Supplementary Data S1, S2; Supplementary Figs. S3, S4). Based on morphological characteristics, bison specimens from Beacon Island (32MN234) and Rustad (32R1775) belonged to B. antiquus 7,24 . The bison remains from other sites were identified as B. bison 25 .
We performed Bayesian phylogenetic and Bayesian skyline analyses with multiple sets of complete mitochondrial genomes. First, we built an uncalibrated phylogenetic tree for a bison global population using available mitogenome datasets, consisting of our 40 mitogenomes from Holocene B. bison and B. antiquus previously radiocarbon dated 19 , 20 mitogenomes of B. priscus of the best quality without stretches of undetermined nucleotides from radiocarbon dated remains originating from eastern and western Beringia 4,27 , and 12 modern B. bison representing the currently existing modern mtDNA haplotypes 28 . The phylogenetic analysis indicated strong statistical support for separating a group consisting of Holocene and modern bison and three haplotypes (KX269120, KX269121, and KX269143) from B. priscus remains excavated in Clover Bar Sand and Gravel site in Alberta and corresponding to Clade 1A in 22 . The sister clade to the above group was composed of two sequences from Alaskan Beringia dated to the terminal Pleistocene. These two mtDNA sequences were isolated from B. priscus specimens from Inuvik, Northwest Territories (KX269125) and Old Crow, Yukon (KX269139) and dated to the time when the corridor between Laurentide and Cordilleran ice sheets opened for migrations after the ice sheets that covered the area 23,000-13,400 ya retreated 17 (Fig. 2A).
Next, we used radiocarbon dates to calibrate the Bayesian phylogenetic tree for a late Pleistocene-Holocenemodern set (Fig. 2B)  , which correlates with the end of the LGM. It also coincides with the temporal range determined as the time of a common ancestor (15,000-23,000 ya) for fossil bison that lived south of 60° N and modern B. bison determined by the mtDNA control region 1,17 . This was a period when the corridor between the Laurentide and Cordilleran ice sheets was closed to migrations 17  Diversification of mtDNA haplotypes. We determined that the most active diversification of mtDNA haplotypes in the NGP occurred 12,000-8000 ya during the early Holocene when bison moved to occupy new territory of the retreated ice sheets and drained Lake Agassiz 29 , as grasslands replaced spruce forests. The calibrated Bayesian tree shows 22 nodes of the splitting events in the early Holocene versus 8 during the late Pleistocene-Holocene transition, 7-in the mid Holocene, and 14-in the late Holocene ( Supplementary Fig. S5). This result demonstrates that Holocene and modern mtDNA lineages originated in different periods, with most haplotypes appearing in the early Holocene (20 haplotypes, 36.4%) and late Holocene (23 haplotypes, 41.8%) (Supplementary Data S3). Many mtDNA haplotypes that originated in the temporal period from the terminal We also confirmed a close relationship between the mtDNA haplotypes of B. antiquus and B. bison indicated previously through analysis of mtDNA control region 6 . This observation is in agreement with the conclusion that B. antiquus did not become extinct but instead evolved to B. bison through phenotypic and morphological adaptation to changing environments 7 . Our data also corroborate a close relationship of Holocene mtDNA lineages from the NGP with mtDNA from modern samples and B. priscus from the northmost part (Alberta) of the Great Plains.
Demographic analysis. Demographic analysis based on the Late Pleistocene-Holocene-modern set of mitogenomes revealed that long-term population stability occurred from 180,000 ya through 25,000 ya with an N e of 24,376 (95% credibility interval of 11,233-48,188) before the decline and the maximum N e of 26,548  4,27 , and 12 bison representing the primary modern mtDNA haplotypes 28 . Different groups of bison samples are labeled by the following colors: purple for the Holocene B. bison (this study), green for B. antiquus from Rustad (this study), blue for B. antiquus from Beacon Island (this study), orange for B. priscus 4,27 , and black for modern bison 28  www.nature.com/scientificreports/ (15,244) ~ 57,000 ya. A sharp population decline began ~ 20,000 ya during the LGM, dropping almost tenfold to N e of 2495 (985-7031) around 10,000-12,000 ya. This period marked the extinction of many megafaunal taxa in North America, including B. priscus in the Arctic. After 10,000 ya, the population began to rapidly grow, reaching the pre-LGM level with N e of 25,709 (12,979) in the mid Holocene (5500 ya). We then detected a slow decline of bison population beginning in the late Holocene approximately 4000 ya through recent times when N e fell to 13,355 (3726-45,151) (Fig. 3A). www.nature.com/scientificreports/ In contrast, demographic analysis of NGP bison in the Holocene based on mitogenome sequences from 40 specimens with the temporal range of 12,226-167 years before present (yr BP) revealed that population dynamics were distinct in several points. This bison population began its origin near the retreating Laurentide ice sheet that covered northern and eastern North Dakota and eastern South Dakota as recently as 12,300 ya 26 (Fig. 1) and collapsed in the nineteenth century due to the anthropogenic bottleneck. We detected three primary dynamic periods for Northern Plains bison. First, population growth from the small N e of 791 (193-3101) began 14,000 ya with melting the ice sheet and the replacement of spruce forests by the grassland after the LGM and the colonization of new territories that became available after the retreating ice sheet and continued through part of the mid Holocene. Then, demographic stability occurred from 6700 to 2700 ya with N e ranging from 16,414 (3009-52,106) to 18,700 (8043-53,302). Finally, populations fell to N e = 8480 (1570-44,171) just before the nineteenth century extermination of the species (Fig. 3B).
We also compared side-by-side the finer population dynamics of the global bison population and the NGP bison population in the last 14,000 years ( Supplementary Fig. S6). This comparison detected a few principal differences. First, despite the decline of the global population that continued until 11,000 yr BP, the Northern Plains bison demonstrated a slow growth since the beginning of colonization of new territories throughout 11,000 ya. In that point the difference in the N e between the global and NGP populations was minimal, 2887 and 1360, respectively. Second, starting ~ 11,000 ya both populations indicated a stable growth till 5500 yr BP with a higher rate in the global population that can be explained by larger territory with various productivity of different ecosystems. The decline first began in the global population shortly after ~ 5500 yr BP while the NGP population was relatively stable until ~ 2700 ya. Around 1800 ya, the difference between both N e s was relatively small with the values of 19,314 and 15,853 in the global and NGP populations, respectively. Both populations were relatively stable before their extermination in the nineteenth century.

Discussion
Our calibrated Bayesian phylogenetic tree consisting of only 40 Holocene mitogenomes confirmed that the common maternal ancestor of Northern Plains mitochondrial genomes lived in the terminal Pleistocene (14,095-23,230 ya) and the main diversification of mtDNA haplotypes occurred in the early Holocene (Supplementary Fig. S7). These haplotypes, without considering indels, indicated haplotype diversity of 0.994 ± 0.008, pairwise difference of 11.25 ± 5.21, and nucleotide diversity of 0.0007 ± 0.0004 (Supplementary Table S1). Some mtDNA sequences from the same site and temporal periods grouped together such as B59 (12,226 yr BP) and B60 (12,109 yr BP) from Beacon Island, B20 (8109 yr BP) and B68 (8114 yr BP) from Rustad. However, all other clades are composed of sequences representing different sites and temporal periods. This analysis disclosed the continuity of diverse mtDNA lineages through the Holocene and the lack of association of mtDNA sequences with archaeological sites in the Northern Plains. The scenario of bison global population dynamics we infer through the Bayesian skyline analysis, based on the largest set of complete mitochondrial genomes analyzed to date, dramatically differs from the demographic model for Beringia bison 1 but is similar in some points to population changes revealed in other studies 10,16 . Differences can be attributed to variation in modeling methods, the addition of new Holocene mtDNA sequences in our study, and a larger number of genetic polymorphisms in complete mitochondrial genomes versus control region sequences. Regardless, we detected the population decline that continued in both models from 25,000 to 23,000 ya until 12,000-10,000 ya achieving a similar N e . After the bottleneck at the terminal Pleistocene, the next wave of population growth and decline happened, indicating another difference in both our and 16 models. In contrast to 16 , we revealed population growth matching pre-LGM levels in the mid Holocene 5000 ya (Fig. 3A), which could be explained by the better survival of bison and more sustainable environments in the Great Plains than in parts of eastern Beringia, including Interior Alaska and Yukon.
Interestingly, the close values of N e we received for the late Pleistocene and the early-mid Holocene may indicate similar vegetative productivity for Pleistocene steppe-tundra ecosystem in Beringia and the early-mid Holocene Great Plains that were able to support large bison populations. Although, the direct comparison of the productivity of both ecosystems is unknown to us, modern analogues of the Pleistocene steppe-tundra have been identified in the Altai Mountains in southern Siberia and true steppes in Central Asia where other large ungulates persist today 31,32 .
Here, we uncovered that the demographic dynamics of the Northern Plains bison was distinct from both the Beringia population 1,10,16 and the late Pleistocene-Holocene-modern global population (this study). One distinction was a two-fold population increase in the NGP immediately following the LGM in contrast to continued population decline in both the Beringian population and the global population. Three parts of the population curve demonstrated a strong correlation with environmental and anthropogenic factors that acted in the Northern Plains after the LGM. The original period of population growth correlated with bison colonization of new lands released by melting ice sheets and the development of grasslands at the beginning of the Holocene (12,000 ya). This original population growth accompanied by the adaptation of bison to climatic and environmental changes is thought to be expressed in the decline in bison body size and mass from 12,500 to 9250 ya 7 . The initial entrance of Paleo-Indians into the territory of present North Dakota 11,500 ya 33 or 9500 ya 34 did not reverse or delay bison population growth.
The bison population in the Northern Plains in the mid and late Holocene had been stable for at least 4000 years (6700-2700 ya) despite climate fluctuations leading to periodic droughts. The longest drought period known as the Altithermal dry and warm period occurred 8000/7000-5000/4000 ya 33,34 . The stability of bison population during this period can be explained by ecological plasticity and adaptation to heterogeneous landscapes in warmer and drier local climate as well as to the human de-population of the NGP 19,33,35 . In contrast to popular opinion, our analysis indicates that American bison were not at their peak in the Great Plains when www.nature.com/scientificreports/ the Europeans first arrived in the sixteenth century 14,36 . Instead, we revealed relatively small but continuous population decline starting 2700 ya. The beginning of mass harvesting of bison 2800 ya based on the application of fire and new weapons and hunting technologies in the late Plains Archaic cultural period 13,33,34,37 is the most parsimonious reason correlating with this decline in the late Holocene. Although mass hunting affected the abundance of bison, it did not result in extirpation, as experienced later in the nineteenth century due to the colonization and industrialization of the Midwest. Our study has revealed considerable diversity of mtDNA lineages during the Holocene. Despite close phylogenetic association, we found no historic mitochondrial genomes matching modern bison 28 . The scope of lost genetic diversity of the Holocene bison during the anthropogenic bottleneck remains unclear as well as the proportion of the survived Holocene lineages. However, considering the complete extermination of bison in the NGP, 100% of these mtDNA lineages were erased in this geographic population and not found in the modern dataset 28 .
Our research has highlighted the importance of historic and geographical populations as indicators of bison community ecology and adaptation to changing environments. Genetic and isotopic evaluation of diverse temporal and spatial bison populations that occupied North America through the Holocene would benefit our understanding of the species in a natural context and enable development of comprehensive strategies to restore bison herds in the future. To best enable bison success, efforts should focus on identifying prairie territory having productivity comparable to the early-mid Holocene Plains and mammoth steppe-tundra.  19 . We cut or drilled a section of teeth or bones for ancient DNA analysis before sampling tissues for radiocarbon dating and stable isotope profiling 19 . Information about archaeological sites, types of specimens, and calibrated age in yr BP was included in Supplementary Data S1. Examples of bison samples are shown in Supplementary Figs. S1 and S2.

DNA isolation.
The project was carried out in a laboratory designated for ancient DNA work at the University of North Dakota. All samples were cleaned by removing soil and dirt deposits using sandpaper. Specimen surfaces were de-contaminated by consequently washing in or wiping with a bleach solution (2% sodium hypochlorite), molecular biology grade water and ethanol. Finally, samples were dried under UV irradiation. All de-contamination steps were conducted within a biosafety cabinet assigned for processing ancient bison materials. We sampled bone and dentin powder samples using a drill with hole saw bits at low drill speeds as recommended elsewhere 38 . We performed DNA isolation from 268 to 370 mg of bone or tooth powders using a DNA extraction method optimized for low-quality samples 39 . Concentrations of DNA isolated from samples were measured by Qubit fluorometry using a high sensitivity assay for low DNA concentrations (Fisher Scientific) (Supplementary Data S1).
Library preparation. We prepared dual-indexed DNA libraries from single-stranded DNA (ssDNA) using the Accel-NGS Methyl Seq DNA library kit with an uracil-tolerant polymerase (Swift Biosciences) according to the manufacturer's protocol for the retention of small DNA fragments. The bisulfite conversion step included in the kit was not used to obtain DNA libraries. The Accel-NGS Methyl Seq DNA library technology adds a low complexity polynucleotide tail with an average length of 8 bases to the 3ʹ end of each fragment and may also add small tails to the 5ʹ ends during end repair steps (Swift Biosciences Technical Note). Depending on the input amount of DNA, the variable numbers of PCR cycles were used to generate DNA libraries according to the manufacturer's protocol. Retention of small DNA fragments was regulated by different ratios of magnetic SPRI beads (Beckman Coulter) in solutions during purification steps, as recommended by Beckman Coulter and Swift Biosciences. Concentrations of DNA libraries were measured by Qubit fluorometry using a high sensitivity assay for low DNA concentrations (Fisher Scientific) (Supplementary Data S1).
The DNA libraries were captured by Daicel Arbor Biosciences (https:// arbor biosci. com/) using a myBaits Expert Mito kit targeting the full B. bison mitochondrial genome (NC_012346). One round of hybridization capture was performed using 4-plex reactions following the myBaits manual (version 4.01) with hybridization and washes steps at 60 °C. Following capture, the first 24 samples were amplified for 10 cycles and the second 23 samples were amplified for 14 cycles. The DNA libraries were sequenced using the Illumina MiSeq platform with SE 150 (150 bp single-end) chemistry. The DNA libraries were sequenced at the University of North Dakota's Genomics Core (24 libraries) and Daicel Arbor Biosciences (23 libraries).

Estimation of consensus mitochondrial genome sequences and summary population statistics.
To produce consensus sequences of the mitochondrial genomes, we tested two approaches. First, DNA reads were processed with BBDuck (http:// sourc eforge. net/ proje cts/ bbmap/) incorporated into Geneious Prime 2021.1.1 (https:// www. genei ous. com) with default settings for trimming the Illumina adapters and low complexity sequences added to the 3ʹ and 5ʹ ends, together with low quality trimming and removing short (< 20 bp) and duplicate reads. DNA reads were mapped using Geneious Prime with default parameters for a modern bison mitochondrial genome (GU946990) sequenced from an animal from a private herd in Montana 28 . To produce consensus mitogenome sequences, we used Geneious Prime with consensus threshold set at the highest quality of 95%. www.nature.com/scientificreports/ Second, because BBDuck does not include the option of removing PCR duplicates, we processed the raw DNA reads and the DNA reads mapped to GU946990 using Paleomix 40 . Short reads with length less than 25 bp were discarded with AdapterRemoval. Picard MarkDuplicates tools v2.26 incorporated in this pipeline were set to "filter" to remove PCR duplicates. For consensus SNP calling, we selected the depth of coverage equals 3 40 .
In total, we obtained 40,860,967 raw DNA reads. The result of reads' processing in Geneious and Paleomix is shown in Supplementary Data S1. The estimated mean coverage ranged from 24 to 7775 in Geneious and from 19 to 299 in Paleomix. The total number of reads mapped to GU946990 ranged from 2818 to 930,448 in Geneious and from 2434 to 32,734 in Paleomix. More robust data received in Paleomix were used for final consensus base calling.
To identify single nucleotide polymorphisms (SNPs), the consensus mitochondrial genome sequences obtained in this study were aligned with the reference bison mitogenome (GU946990) using MEGA7 41 Table S1).
The estimation of DNA damage parameters. Ancient DNA damage parameters of single-end DNA reads were analyzed by mapDamage 2.0 43 . In ancient DNA molecules, it is expected that single-end DNA reads have an increased frequency of C to T misincorporations at the 5ʹ ends as the consequence of cytosine deamination ( Supplementary Fig. S3) and higher frequency of guanine or adenine in the − 1 position of DNA breaking points in front of the 5ʹ end of DNA reads (Supplementary Fig. S4) 44 .

Phylogenetic analyses.
To understand the evolutionary relationships of 40 Holocene bison mitogenome sequences obtained in this study, we selected the best 19 mitochondrial genome sequences of B. priscus from 4 and one-from 27 isolated from dated remains and having no stretches of undetermined (N) nucleotides. We also included 12 modern bison mitogenome sequences representing the main modern haplogroups 28 .
We have excluded from consideration partial mitogenome sequences with missing parts from Holocene steppe bison from the Yukon 21 and from two historic bison remains sampled in northern Wyoming in 1856 and southern Montana in 1886 45 .
Bayesian phylogenetic analysis of mtDNA sequences with MrBayes. We utilized MrBayes 3.2. 6 46 implemented in Geneious Prime to estimate a Bayesian phylogenetic tree using the HKY + I + G substitutional model as determined by jModelTest2 47 . We ran the Markov chain Monte-Carlo (MCMC) analysis using 1,100,000 generations with 4 heated chains, a burn-in length of 100,000, and sampling every 200 generations. The MCMC convergence was diagnosed by the average standard deviation of split frequencies (ASDSF) 46  Bayesian phylogenetic analysis of mtDNA sequences with the median calibrated radiocarbon dates used to calibrate the tree. We analyzed two datasets of bison mtDNA sequences. The first set included complete mitochondrial genomes from 40 Holocene bison (this study), 20 Pleistocene B. priscus (19 sequences from 4 and one sequence from 27 ) and 12 modern bison lineages 28 . The second set contained only 40 Holocene bison mitogenome sequences obtained in this study.
We used BEAST 1.10.4 49 to build Bayesian phylogenetic trees based on complete mitochondrial DNA sequences and their radiocarbon age using the HKY + I + G substitution model for the late Pleistocene-Holocene-modern set and the HKY + G substitution model for the Holocene set according to the jModelTest 2 definition 47 , a strict molecular clock, and the Bayesian Skyline coalescence model with ten groups. The MCMC analysis ran 10,000,000 generations with sampling each 1000 generations with the convergence diagnosed in Tracer 48 .
TreeAnnotator (http:// beast. bio. ed. ac. uk) was used to make calibrated phylogenetic trees. TreeAnnotator in BEAST 1.10.4 package 49 was used to compile the sampled trees in a single target tree to annotate summarized information. Phylogenetic trees were visualized using FigTree 1.4.4 (http:// tree. bio. ed. ac. uk/ softw are/ figtr ee/). Demographic analysis. We analyzed two datasets of bison mtDNA sequences described above using the Bayesian Skyline model previously applied to steppe bison demographic studies 16,50,51 . To reveal bison population demography over time, we analyzed log files generated in BEAST 1.10.4 using Tracer 1.7.1 48 . Tracer produced graphical plots of demographic reconstructions showing changes of N e τ (the y axis), where N e is an effective population size based on mtDNA corresponding to female effective population size and τ is a generation time, over time in radiocarbon years before the present (the x axis) 16 . We selected 8 years as the bison generation time in this study 30 although bison generation time has ranged 3-10 years by different authors 7 .

Data availability
All data are available in the main text or the Supplementary materials. Original fastq files with accession number PRJEB62763 are deposited in the European Nucleotide Archive (ENA).